#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use db_parameters;
use ensembl_parameters;
use chromosomes;

use DBI;


print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_prot_exons

	Description:
		fills
		taxon.prot_exons
		taxon.gene_loci

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE


#_________________________________________________________________________________________

#	retrieve_chromosome_lengths

#_________________________________________________________________________________________
sub retrieve_chromosome_lengths($\%)
{
	my ($dbh, $chrm_lengths) = @_;
	print STDERR "\tReading chromosome lengths...\n";
	my $data = $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
        species,
		name,
		length
	FROM
		taxon.chromosomes
PL/SQLCMD

	for my $row(@$data)
	{
		$chrm_lengths->{$row->[0]}{$row->[1]} = $row->[2];
	}
}	

#_________________________________________________________________________________________

#	retrieve_prot_id_to_gene_id

#		Do this here rather than ens_id to avoid 
#		order dependencies and minimize coupling

#_________________________________________________________________________________________
sub retrieve_prot_id_to_gene_id($\%)
{
	my ($sequence_dir, $prot_id_to_gene_id) = @_;
	
	print STDERR "\tReading prot id to gene identifier\n";
	my $ens_id_file = "$sequence_dir/prot_id_to_gene_id.map";
	die "Error:\n\t$ens_id_file does not exist\n" unless (-f "$ens_id_file");
	open ENS_ID, $ens_id_file or die;
	my $cnt_line = 0;
	while (<ENS_ID>)
	{
		++$cnt_line;
		my (undef, $gene_id, $prot_id, undef) = split /\t/;
		$prot_id_to_gene_id->{$prot_id} = $gene_id;
	}
}




#
#	connect to panda
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);


#
#	get ensembl taxa
#
my %ensembl_taxa;
for my $ens_data(get_ensembl_taxon_data())
{
	++$ensembl_taxa{$ens_data->[TAXON]};
}

#
#	Get chromosome lengths
#
my %chrm_lengths;
retrieve_chromosome_lengths($dbh, %chrm_lengths);


#
#	get custom
#
my %taxon_data = get_taxon_data;
for my $taxon(keys %taxon_data)
{
	next if exists $ensembl_taxa{$taxon};

	print STDERR "\tProcessing non-Ensembl species $taxon...\n";
	my $sequence_dir = $taxon_data{$taxon}[SEQ_DIR];


	#
	#	Get prot id to gene id lookup. Do this here rather than ens_id to avoid 
	#		order dependencies and minimize coupling
	#
	my %prot_id_to_gene_id;
	retrieve_prot_id_to_gene_id($sequence_dir, %prot_id_to_gene_id);

		
	my %gene_id_loci;
	
	
	#
	#	prot_exons
	#
	# open file
	use constant CHROMOSOME => 0;
	use constant STRAND		=> 1;
	use constant START		=> 2;
	use constant FINISH		=> 3;
	use constant GENE_ID	=> 4;
	
	print STDERR "\tReading exon loci\n";
	my $loci_file = "$sequence_dir/exon.loci";
	die "Error:\n\t$loci_file does not exist\n" unless (-f $loci_file);
	open LOCI, $loci_file or die;
	my $loci_data_convention = <LOCI>;
	chomp $loci_data_convention;
	my $andreas_mode = $loci_data_convention eq 'ANDREAS' ? 1 : 0;
	print STDERR "\tUsing Andreas data (backwards for -ve strand)\n" if $andreas_mode;
	
	my $id_prefix = $taxon_data{$taxon}[ID_PREFIX];
	print STDERR "\tClear old data matching '$id_prefix'\n";
	$dbh->do("SET WORK_MEM=1000000;");
	$dbh->do("DELETE FROM taxon.prot_exons ".
			 "WHERE prot_id ~ '^$id_prefix'");
	$dbh->do("COPY taxon.prot_exons".
				"(prot_id, exon_id, start, finish) FROM STDIN");

	my $cnt_exons = 0;
	my $cnt_short_exons = 0;
	
	while (<LOCI>)
	{
		chomp;
		my ($prot_id, $chr, $strand, $exon_id, $beg, $end) = split /\t/;
		die $_ unless defined $end;
		if ($andreas_mode)
		{
			die $_ unless exists $chrm_lengths{$taxon}{$chr};


			# Ensembl coordinates should be inclusive beg < end from 1
			# Andreas coordinates non-inclusive from 0
			# Also counts backwards if -ve strand
			#
			# len = 10
			#
			# andreas 3-5 +strand
			# 
			# 0123456789
			#    **
			# 1234567890
			# Ensembl 4-5
			#
			#	Begin +=1, end stays the same
			#
			#
			# andreas 3-6  -strand
			#
			# 9876543210
			#     ***
			# 1234567890
			# Ensembl 5-7
			#
			# beg = len - end + 1	10 - 3       = 7
			# end = len - beg       10 - 6 + 1   = 5
			if ($strand)
			{
				++$beg;
			}
			else
			{
				($beg, $end) = ($chrm_lengths{$taxon}{$chr} - $end - 1, 
								$chrm_lengths{$taxon}{$chr} - $beg);
			}
		}
		#print "$_\n" if $beg == $end;
		++$cnt_short_exons if $end - $beg + 1 <= 10;
		++$cnt_exons;
		$dbh->pg_putline("$prot_id\t$exon_id\t$beg\t$end\n");

		next unless exists $prot_id_to_gene_id{$prot_id};
		my $gene_id = $prot_id_to_gene_id{$prot_id};
		

		#
		#	gene loci is the greatest span over all the exons
		#
		$strand = $strand ? 't' :'f';
		if (exists $gene_id_loci{$chr}{$gene_id})
		{
			die unless $gene_id_loci{$chr}{$gene_id}[CHROMOSOME] eq $chr;
			die unless $gene_id_loci{$chr}{$gene_id}[STRAND] eq $strand;
			if ($gene_id_loci{$chr}{$gene_id}[START] > $beg)
			{                      
				$gene_id_loci{$chr}{$gene_id}[START] = $beg;
			}                      
			if ($gene_id_loci{$chr}{$gene_id}[FINISH] > $end)
			{                      
				$gene_id_loci{$chr}{$gene_id}[FINISH] = $end;
			}
		}
		else
		{
			$gene_id_loci{$chr}{$gene_id} = [$chr, $strand, $beg, $end];
		}
	}
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy;

	print STDERR "\t\t$cnt_exons exons entered\n";
	print STDERR "\t\t$cnt_short_exons exons entered < 10 bp\n";
	

	
	#
	#	gene_loci
	#
	print STDERR "\tSort gene loci by chromosomes and order\n";
	my @chromosomes = sort_chromosomes(keys %gene_id_loci);
	# put into list and sort by location
	my @gene_loci;
	for my $chrm(@chromosomes)
	{
		print STDERR "\tchromosome $chrm:\t";
		my @chrm_gene_loci;
		for my $gene_id(keys %{$gene_id_loci{$chrm}})
		{
			push(@chrm_gene_loci, [@{$gene_id_loci{$chrm}{$gene_id}}, $gene_id]);
		}
		@chrm_gene_loci = sort {($a->[START] + $a->[FINISH]) <=> ($b->[START] + $b->[FINISH])} @chrm_gene_loci;
		print STDERR scalar @chrm_gene_loci, " genes\n";;
		push(@gene_loci, @chrm_gene_loci);
	}
	print STDERR "\t", scalar @gene_loci, " genes in total\n";;
	my $db_origin = $taxon_data{$taxon}[PANDA_DB_ORIGIN];
	
	
	print STDERR "\tRemoving old data from panda\n";
	$dbh->do("DELETE FROM taxon.gene_loci ".
			 "WHERE db_origin = '$db_origin'");
	print STDERR "\tCopying gene loci to panda\n";
	$dbh->do("COPY taxon.gene_loci".
				"(chromosome, strand, start, finish, gene_id, ".
				"db_origin, gene_position) FROM STDIN");
	my $gene_pos = 0;
	for my $locus(@gene_loci)
	{
		$dbh->pg_putline(join ("\t",  @$locus, $db_origin, ++$gene_pos)."\n");
	}
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy;
}
print STDERR "\tVacuum analyse taxon.prot_exons\n";
$dbh->do("SET MAINTENANCE_WORK_MEM=500000;");
$dbh->do("VACUUM FULL taxon.prot_exons;");
$dbh->do("VACUUM FULL taxon.gene_loci;");

print STDERR "\tCompleted\n";

log_pipeline_stage($dbh, 'NULL', 20, 'table prot_exons', $start_time);
$dbh->disconnect();

